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■ We reanalyze the results of our extensive Hartree-Fock + BCS calculation from new 

' points of view paying attention to the properties of unstable nuclei. The calculation has 

^vQ , been done with the Skyrme SIII force for the ground and shape isomeric states of 1029 

even-even nuclei ranging 2 < Z < 114. We also discuss the advantages of the employed 
three-dimensional Cartesian-mesh representation, especially on its remarkably high preci- 
sion with apparently coarse meshes when applied to atomic nuclei. In Appendices we give 
the coefficients of finite-point numerical differentiation and integration formulae suitable for 
Cartesian mesh representation and elucidate the features of each formula and the differences 
from a method based on the Fourier transformation. 
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! SI. Introduction 

oo 

■ Among theoretical attempts which aim at treating all the nuclides in a single 

i framework, the simplest category seems to be the mass formula. Hence, let us focus 

^ I on the nuclear masses in the first place, although the scope of this paper is not 

• restricted to the masses. The most familiar mass formula is certainly that of Bethe 

and Weizsacker,^''^ which expresses the nuclear masses as a function of the number 
of neutrons (A'') and that of protons (Z) as, 

^' E{N,Z) = ayA + asA^/^ + ai{N - ZfA-^ + acZ^A-'^/'\ (1-1) 

where A = N + Z is the mass number. The four terms on the right-hand side 
of this equation are called respectively the volume, the surface, the symmetry and 
5^ ■ the Coulomb terms. These terms come from the liquid-drop picture of the atomic 

^ nucleus, in which the density distribution of nucleons is assumed to be spatially 

uniform inside nuclei and the uniform density value is common to all the nuclei. 

Incidentally, it is worth stressing that even the form of the symmetry term is 
completely understandable within the liquid-drop picture. It is not necessary to 
bring up the Fermi gap model like, e.g., in p. 5 of Ref.^' The form of this term 
can be naturally derived for a liquid drop by admitting that the coefficient of the 
volume term should have a dependence on the ratio of the constituents. Because one 
knows that the volume energy becomes maximum when the proportion of protons 
(or neutrons) is 50%, the volume term should be modified as, 

avA («v + aVQ-f)'U = ay A + ^-^l^^^^ , (1-2) 
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if one retains terms only up to the second order in Z/A. On the right-hand side 
of this equation, the first term is the original volume term and the second term 
is identical to the symmetry term by rewriting = aj. If one applies a similar 
consideration to the surface term one obtains a so-called surface symmetry term, 
a[{N — This term is sometimes added to the liquid-drop mass formula of 

Eq. recently. 

Now, by changing the values of the four coefficients (ay, as, oi and ac), the 
root-mean-square (r.m.s.) error from all the experimental data of about 2,000 even- 
even nuclei available at presenP^'l^ can be decreased down to 3.5 MeV. Inclusion 
of the surface symmetry term can decrease the error further. The accuracy of 3.5 
MeV is very small compared with the binding energies of heavy nuclei (~ 1 GeV). 
However, the requirement for theoretical predictions of the energies of unknown 
nuclei is something more precise. The separation energy of a neutron, S^, decreases 
as the neutron number is increased. Its changing rate can be roughly estimated 
to be 

of which the right-hand side of the first (approximate) equality is the reciprocal of 
the single-particle level density at the Fermi level obtained by assuming a harmonic 
oscillator potential with ^osc = 41^~^/^ MeV. The right-hand side of the last 
equality, which is the expression for nuclei on the neutron drip line, becomes as small 
as 100 keV for the heaviest nuclei. It means that the precision of mass predictions 
must be 100 keV in order to predict the location of the neutron drip line. 

The first step to decrease the error is to take into account the shell effect. For 
example, the TUYY mass formula-' achieved an r.m.s. mass error of 538 keV. The 
number of fitting parameters are 6 for the gross part (corresponding to the parame- 
ters of the Bethe-Weizsacker formula) while that for the shell part is as many as 269: 
There is one parameter for each value of Z in an interval 1 < Z < 112 and one for 
each value of in 1 < A'^ < 157. Generally speaking, less number of parameters are 
preferable for the reliability of the extrapolation to nuclei not synthesized yet. One 
usually switches to less phenomenological models in order to reduce the number of 
parameters. 

What should be considered next is that the shell effect is strongly dependent 
on deformation. Indeed, the research group which presented the TUYY mass for- 
mula has started the development of a new method which includes the concept of 
deformation.-^'"^' Of course there are different directions to proceed. For example, a 
model based on the shell-model configuration mixing has achieved an r.m.s. error of 
375 keV with 28 parameters.!^ The latest good review on the methods of theoretical 
predictions of nuclear masses can be found in Ref.^ 

The most elaborately developed model which takes into account deformation 
seems to be the finite-range droplet model with a microscopic shell correction (FRDM), 
whose latest update was done by Moller et al.Bl'' Another extensive calculation 
was carried out by Aboussir et al.^ in the extended Thomas-Fermi plus Strutinsky 
integral method (ETFSI). 

These two methods can be regarded as approximations to the Hartree-Fock (HF) 
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method. The straight-forward solutions of the HF equations including deformation 
require long computation time even with present computers in order to cover thou- 
sands of nuclei. We understand that the first calculation of such a large scale which 
have been published is our calculatior^^ done in the framework of the HF-I-BCS 
method with the Skyrme SIII force. After the publication of our paperj^^l only a 
few similar extensive calculations have been done. Tondeur et al. have performed a 
set of such extensive calculations in the framework of the HF-I-BCS and proposed 
a new Skyrme force parameter set MSk7 best suited for their truncated oscillator 
basis.^'^ Another attempt have been carried out in the relativistic mean-field 
framework.^ 

In section 2, we discuss about the Cartesian mesh representation of the wave 
functions, which is a feature of our calculations. Details not included in Ref."^ are 
relegated to the Apppendices. 

In section 3, we reanalyze the results of our extensive calculation in order to give 
discussions from new points of view and to compare the results with those of other 
models. 

§2. Cartesian-Mesh Representation 

The most usual method to solve the mean-field equations without assuming the 
spherical symmetry is to expand single-particle wave functions in terms of truncated 
anisotropic harmonic-oscillator basis.^'^ Instead, we employ a three-dimensional 
Cartesian-mesh representation^^ We put a nucleus in a box containing ~ lO'^ mesh 
points and express each single-particle wave function in terms of its values at the 
mesh points. The advantages of this representation can be summarized as follows: 

(1) It has no prejudice concerning the shape of the nucleus. On the other hand, 
in methods based on expansions in the anisotropic oscillator basis, the shape of 
the solution should be similar to that of the basis in order that the truncation 
of the basis does not affect the solution. To fulfill this requirement, one usually 
solves the mean-field equation imposing constraints on the quadrupole moments 
such that the resulting shape agrees with the anisotropy of the basis. Then, to 
obtain the ground state, one has to optimize the anisotropy of the basis such that 
the energy of the solution is minimized. This procedure is not only cumbersome but 
also unmanageable in treating exotic states: For example, if the protons and the 
neutrons may have different shapes, the dimension of the parameter space for the 
optimization is squared. On the other hand, in the mesh representation, one can 
treat various shapes with the same mesh. Consequently, the optimization procedure 
is not necessary. 

(2) In the mesh representation the asymptotic form of the wave function at large 
radius can have arbitrary form, while in the oscillator basis representation it must 
have a Gaussian tail and therefore loosely bound single-particle wave functions like 
neutron halcP cannot be correctly described. 

(3) Among systems of many particles found in nature, the atomic nucleus is 
a very suitable object to apply the Cartesian-mesh representation: Inside nuclei, 
the density is roughly constant and therefore the local Fermi momentum is also 
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constant everywhere. This situation favors an expansion in the plane wave basis, 
to which the Cartesian mesh representation can be an accurate approximation by 
using appropriate approximation formulae for derivatives and integrals as discussed 
in Appendices. Consequently, apparently coarse meshes can result in unexpectedly 
high precision: With mesh spacing a ~ 1 fm, there are only a few mesh points 
in the surface region. Nevertheless, the relative error of the total energy and the 
quadrupole moment with this mesh spacing is as small as 0.5 % for heavy nuclei. 
The left-hand portion of Fig. ^ shows how the total energy of ^^'^Er changes as a 
function of the mesh spacing. At a=l fm, the error is ~7 MeV, which is indeed small 
(0.5%). Furthermore, as shown in the right-hand portion of the figure, the energy 
difference between prolate and oblate solutions is by far more precisely determined 
(error ~ 0.4 MeV) with a=l fm. Owing to this independence of the error from 
the shape, we have found it possible to correct the total energy by adding a simple 
function of only Z and N to decrease the error down to ~0.2 MeV 
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Fig. 1. Left-hand portion shows the total energy of ^'^"Er (for the prolate minimum, the oblate 
minimum, and the spherical saddle point) versus mesh spacing a. The circles, triangles, and 
squares are the results of calculations while the curves represent functions E = cq + c-^c^ — c^a^ 
where a are determined through least-square fittings. Right-hand portion gives the energy 
difference between prolate and oblate minima for five nuclei. 



(4) The algorithm of the calculation is simple. This simpleness increases the 
reliability of computer programs based on the algorithm. It also makes the programs 
efficiently executable with vector-processor type computers. The programs can also 
be fitted to each unit of massively parallel computers because the necessary memory 
is not huge. 



§3. HF-I-BCS Calculations with the Skyrme Force 



Using the Cartesian mesh representation, we have calculated the ground and 
shape isomeric states of even-even nuclei with 2 < Z < 114 and A'^ ranging from 
outside the proton drip line to beyond the experimental frontier in the neutron-rich 
side. We have obtained spatially localized solutions for 1029 nuclei, together with 
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the second minima for 758 nuclei. Details not described in this paper can be found 

The Skyrme is an effective interaction widely used in mean-field cal- 

culations. It is essentially a zero-range force but modified with the lowest order 
momentum dependences to simulate the finite-range effects, a density dependence to 
reproduce the saturation, and a spin-orbit coupling term. The is one of the 

many parameter sets proposed for the Skyrme force. It features good single-particle 
spectra and accurate N — Z dependence of the binding energy.'^ 

We have used a computer program named EVSW In the program, one places 
an octant of a nucleus in a corner of a box (13 x 13 x 14 fm^ for Z < 82 and 14 x 
14 X 15 fm^ for Z > 82), imposing a symmetry with respect to reflections in x-y, 
y-z, and z-x planes (the point group D2}^). The mesh spacing is 1 fm as explained 
in the last section. 

For the pairing, we employ a seniority force, whose pair-scattering matrix ele- 
ments are defined as a constant G-r (r distinguishes between proton and neutron) 
multiplied by a cut-off factor which is a function of the single-particle energy e. The 
cut-off is set at e = (Fermi level -|- 5 MeV) with smearing width of 0.5 MeV. For 
neutrons, the cut-off function is multiplied furthermore by 6[—e). The strength Gr 
is determined for each nucleus such that the continuous spectrum approximation us- 
ing the Thomas-Fermi single-particle level density reproduces the classical empirical 
formula Z\ = 12 MeV/\/l. 

Since the pairing correlation has strong influences on deformations, one must 
not trifle the choice of the pairing force if one wants to perform deformed mean-field 
calculations. Our treatment is simple in the sense that it employs the seniority force 
but it is an advanced treatment in the sense that the force strength is determined as 
a function of the size of the configuration space for the pairing correlation. 
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Fig. 2. Dependence of the seniority pairing force strength G renormahzed in such a way that the 
resulting pairing gap agrees with the experimental value (top portion) and the dependence of the 
energy gain due to the pairing correlation calculated with the strength G (bottom portion) on 
the cut-off energy. Single-particle levels of the Nilsson potential are used in the BCS calculation. 

It is important not to include positive energy neutron orbitals in BCS calcu- 
lations. Otherwise the resulting density distribution has unphysical neutron gas 
spreading over the normalization box. For neutron-rich nuclei in which the Fermi 
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level is not far from zero energy, the pairing configuration space restricted to negative- 
energy HF orbitals becomes too small to obtain a reliable result within the BCS 
scheme. 

Figure HI shows the results of a test calculation to illustrate the inadequacy of 
too small pairing spaces for BCS calculations. The calculation is done using the 
single-particle spectrum of the Nilsson model, which does not include the continuum 
spectrum and free from the problem of the formation of the unphysical neutron gas. 
The strength G of the seniority pairing force is determined in such a way that, for 
each value of the cut-off (measured from the Fermi level), the neutron's pairing gap 
Z\n calculated within the BCS scheme agrees with the experimental value of 1.45 
MeV. The value of thus defined strength G is shown in the upper half portion of 
Fig. [21 The bottom half portion shows the energy gain due to the pairing correlation 
calculated with the value of G shown in the top portion. The energy gain does 
not change sizeably for the cut-off energy greater than ~ 4 MeV. For cut-offs less 
than that value, however, the energy gain is non-negligibly smaller than the values 
calculated with enough large cut-offs. 




100 

NEUTRON NUMBER 

Fig. 3. Even-even nuclei covered by our extensive calculation (open circles). Dots are the uncalcu- 
lated even-even nuclei between the drip lines. Two-proton and two-neutron drip lines are drawn 
with solid lines. Dot lines are those from the results of FRDM.^ 



From the above consideration concerning the effect of the continuum states on 
the neutron's pairing correlation, it is not very meaningful to apply the HF-I-BCS 
scheme to nuclei in which the Fermi level of neutrons is higher than ~ —4 MeV. Thus 
we restrict our calculation within only several neutrons beyond the experimental 
neutron-rich frontier. Figure 01 shows thus selected even-even nuclei. The two-proton 
drip line (the upper one of the solid lines) is deduced from the total energies of our 
HF-I-BCS calculation. Solution of nuclei outside the proton drip line can be obtained 
owing to the Coulomb barrier. The two-neutron drip lines (the lower one of the solid 
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lines) is obtained from the Bethe-Weizsacker type mass formula of Eq. whose 
four coefficients are determined through a least-square fitting to the total energies 
from our HF+BCS calculation. One can see that roughly half of the even-even nuclei 
are included in our calculation. 

In light-mass region, the neutron's Fermi levels of some of the calculated nuclei 
are higher than —4 MeV. For such nuclei, the renormalized pairing force strength 
becomes very large. We cut the strength at 0.6 MeV. As a consequence, our solutions 
for such light-mass nuclei near the neutron drip line do not have pairing correlation. 

The current subject of our research is the calculation of the remaining half of 
the nuclear chart on the neutron rich side by developing a feasible Hartree-Fock- 
Bogoliubov framework which enables one to include the effect of the Hartree-Fock 
continuum (states which are in the continuum part of the Hartree-Fock single-particle 
spectrum) on the pairing correlation. We think it preferable if the wave functions 
are expressed in the coordinate space with a three-dimensional Cartesian mesh in 
order to describe possible large surface diffuseness like thick skins and halos as well 
as arbitrary deformations. A method utilizing the localization of HFB canonical 
basis seems to be very promising .-^"^ 

§4. Analyses of the Results 

In this section we reanalyze the results of our extensive HF+BCS calculations,^ 
by presenting the results from new aspects. Note that all the resulting data of our 
calculations are available on the WEB.^^' In principle one can reproduce the figures 
given in this section using only the data on the WEB. 

4.1. Mass 

Concerning the nuclear masses, the r.m.s. error of our result turns out to be 2.2 
MeV, which is not so good as those of the latest mass formulae (typically 0.5 MeV). 
Note that this was the first published estimation of the r.m.s. mass error of any of 
the Skyrme forces. It is an interesting question how much the error can be decreased 
by improving the parameters of the Skyrme force through an extensive fitting to 
currently available experimental masses.^'^ To achieve the precision of 100 keV, 
improvement of the treatment of the pairing is certainly necessary. It may also be 
necessary to add some new terms to the Skyrme force. Some think that corrections 
for correlation energies should be included. 

Leaving this question which is out of the scope of this paper, let us reanalyze 
our result from different aspects not taken in Ref.^ 

The top portion of Fig. |1] shows the errors of the calculated masses as a function 
of the neutron number A^. The bottom portion plots the same data versus the 
proton number Z. Nuclei oi N = Z are designated with open circles (appearing in 
N, Z < 40) while the others with dots. 

One can see that the N = Z nuclei are always at the highest peaks of the 
isotope and the isotone chains to which they belong. This means that our mean- 
field calculation does not have the mechanism corresponding to the Wigner term of 
the mass formulae. It is now widely understood that the Wigner comes from the 
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Fig. 4. Errors of the nuclear masses calculated in HF+BCS with Skyrme SIII force plotted as a 
function of the neutron number TV (top portion) and the proton number Z (bottom portion). 
Isotopes and isotones are connected with lines. The number printed near each line of isotope or 
isotone chain means the proton or neutron number of the chain. 



T = pairing. On the other hand, in the present HF+BCS approach the pairing 
correlation is considered only between like nucleons, which means a T = 1 pairing. 

The behavior of the errors looks very systematic in region Z,N > 40 compared 
with in lighter-mass region. One can see from the top figm'e that mass difference 
becomes minimum (i.e., the calculated values are more bound than experimental 
values) always at neutron's spherical magics of N = 50, 82, and 126. The situation 
is less clear in the bottom figure, but at least a spherical magic of protons, Z = 50, 
is a minimum of the isotone chains except for the N = 82 chain. This observation 
suggests that the employed force has too large a stiffness against deformation. 

Concerning the possibly too large surface tension of the employed force. Fig. [S] 
plots the errors of the calculated nuclear masses divided by A^^^ (which is propor- 
tional to the error per unit surface area) versus i?surf — Ij where i?surf is the ratio of 
the surface area to that of the spherical shape. The increase of the surface energy 
due to deformation is proportional to -Bgurf — 1- The surface area is calculated using 
liquid-drop shape parameters, whose definition in term of moments is given later in 
the next subsection. The error is negative for half of the spherical nuclei (.Bsurf = 1); 
while for deformed nuclei it is roughly constant and independent of .Bsurf (i-e., the 
size of deformation). The situation is also independent of whether the nucleus is light 
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Fig. 5. Errors of nuclear masses divided by a mass-number dependence factor A^^^ as a function 
of the surface factor Bsurf subtracted with 1. Open circles, triangles, and squares are used for 
the sake of distinction of the heaviness of the nuclei. The symbols in the left-hand side of the 
vertical dot line Bsm-f — 1=0 should be regarded as on the line: They are shifted to the left 
only to avoid clustering of too many symbols. 



or heavy as one can see three kmds of symbols (circles, triangles, and squares) mixed 
up in the plot. These results do not suggest that the tendency that deformed nuclei 
have too large masses can be ascribed solely to the surface energy of the employed 
force. 

4.2. Deformation 

Fig. El shows a comparison between the experimental and the calculated intrinsic 
quadrupole moments. For each nucleus whose B(E2; 0"*" — > 2"*") is experimentally 
known (289 even-even nuclei ranging over 4 < Z < 98),'^'^' a dot is marked at a point 
whose abscissa is equal to the experimental value and its ordinate to the theoretical 
value calculated in this paper. The agreement between the experiment and the theory 
is excellent in most cases. Only a few nuclei exhibit large discrepancies. Three such 
nuclei whose element names are indicated in the graph (^'^'^Pt, ^^^Ra, and 222Th) 
are around the shape transitional point of their isotope and isotone chains. It is 
relatively difficult to predict the deformation of shape transitional nuclei. 

It has been reported that our results have roughly the same quality as the 
other principal theoretical approaches concerning the reproduction of experimental 
quadrupole moments.!^ 

Let us define the deformation parameters in terms of multipole moments for 
the sake of comparison with the results of other models like FRDM and also for the 
easiness to imagine the shapes. We define the deformation parameters of a HF-I-BCS 
solution as those of a sharp-surface uniform-density liquid drop which has the same 
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mass moments as the solution has. The mass density of the hquid drop is expressed 



as 



oo I 



p(r) = po 9{R{r) -\r\), R{r) =Ro{l + Y.Yl ^i^^Ur 



(4-1) 



1=0 m=-l 



where 6 is the step (Heaviside) function. The necessary and sufficient conditions 
on aim to fulfill the reality of R{r) and the symmetry are that / and m are 



even numbers and a^^ 



O'l-m- We set aim = for / > 6 and determine the 



remaining seven parameters po, Ro, 0,20, a22, 04O) 0,42-, and 044 such that the liquid 
drop has the same particle number, mean-square mass radius, and mass quadrupole 
{r'^Y2m) and hexadecapole (r'^Y/im) moments as the HF+BCS solution has. 

By using the mass number A and the mass moments, one obtains the defor- 
mation parameters to express the shape of the mass distribution. The deformation 
parameters of proton (neutron) distribution can be calculated in a similar way by 
using proton (neutron) number Z (N) and the moments of proton (neutron) distri- 
bution. 

Let us check the adequacy of our definition of the shape parameters by comparing 
them with model-independent quantities. 

In the left-hand portion of Fig. [7J we show the correspondence between two 
quantities related to the axially symmetric quadrupole deformation. They are 020 
and 6, where the latter is defined as 
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The number of plotted points is 5361. Each point is for the deformation of either 
mass, proton, or neutron distributions of a HF-I-BCS solution among 1029 ground 
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Fig. 7. The correspondence between liquid-drop-model shape parameters and the related model 
independent quantities. The left-hand portion is for the axially symmetric quadrupole deforma- 
tion while the right-hand portion is for the radius. 



and 758 first-excited states. The solid line represents a2o = (^l'^)"'^''^'^! which is the 
leading-order expression for sharp-surface density distributions. This is a satisfactory 
result in the sense that the values of 020 are never very distant from the widely 
used model independent deformation parameter 5. We have least-square fitted a 
polynomial through order three with the leading order coefficient frozen at the above 
value. The resulting polynomial is 020 = —0.47(5^ -1-0.785^. The maximum 

and the r.m.s. deviations from this polynomial are 0.05 and 0.007, respectively. 

In the right-hand portion of Fig. [3 we plot the liquid-drop-model radius Rq 
versus the r.m.s. radius rrms = (r^)-*^/^. The diagonal solid line is Rq = (|)^''^rrms, 
which is a relation expected to a spherical liquid drop. It is again a desirable result 
that the deviations from this line are never be very large. Owing to deformations 
and the surface diffuseness, Rq is likely to be slightly smaller than (-o) ' rj-ms- 

For 

the plotted 5361 points, the maximum and the r.m.s. deviations are 0.3 fm and 0.06 
fm, respectively. 

In Fig. El we compare the axial quadrupole deformation parameter 020 between 
present HF-I-BCS calculation and the results of the FRDM.E* 

A systematic difference which is most easily noticed is found around Z ^ N ^ 40, 
where the present HF-I-BCS calculation tends to predict smaller deformations than 
the FRDM. This is due to the shape coexistence prevailing in this region. 

An example of the shape coexistence is given in Fig. |UJ For ^^Zr, the HF-I-BCS 
method with the SHI force (solid curve) predicts three minima which are energetically 
competing within 0.6 MeV. The order of the energies of these minima can be altered 
easily by changing the force parameter to SkM"'^ (dash curve) or by decreasing the 
pairing gap by 25% (dot curve). 

Another systematic difference which one may notice in Fig. |H1 is in a long and 
narrow region close to the proton drip line with 94 < Z < 102, where the FRDM 
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Fig. 8. Comparison of the calculated quadrupole deformation parameter 020 between HF+BCS 
with Skyrme SIII forc(P and the finite-range droplet model^ The open (solid) circles desig- 
nate prolate (oblate) nuclei, while the diameter of the circles is proportional to the magnitude 
of the deformation parameter. The two-proton and two-neutron drip lines predicted by each 
model are also drawn. 



en 



\ ■ 


1 ' 1 ' // 
/ / 
/ / 

1 r 

■ '/ 




// 

• / / - 


80^ \ 

Zr \ / — 




- , 1 


1 , 1 , " 



0.4 



-0.2 0.2 

QUADRUPOLE DEFORMATION 6 

Fig. 9. The Potential energy curve of *°Zr along the axially symmetric quadrupole deformation 
path. The solid curve is calculated with the SIII force, the dot curve with the SIII force with 
25% smaller pairing gap, and the dash curve with the SkM* force. 



Hartree-Fock+BCS approach to unstable nuclei with the Skyrme force 13 
predicts oblate shapes while our calculations give small-size prolate shapes. 
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Fig. 10. The potential energy curves for Ip shell even-even nuclei from the HF+BCS with the SIII 
force. The energy is measured from that of the spherical solution. For ^^C the results with the 
Skyrme SII force is also drawn with dash line. 



Figure ITUl shows the potential energy curves for 4 < Z, < 6 nuclei, which are 
obtained with a constraint on the axially symmetric mass quadrupole moment Qz- 
The abscissa represents the deformation parameter 5. 

Among our calculations for the 1029 even-even nuclei, ®Be has the largest defor- 
mation (a2o=0.62), which can be explained in terms of the two-alpha-cluster picture. 
In the bottom-left portion of Fig. ^1 one finds the potential energy curve for this 
nucleus, which exhibits that the deep prolate minimum is the only solution (except 
for the extremely shallow one at an oblate shape). 

For ^-^C, the experimental B(E2)| is very large (corresponding to |a2o|=0.59) and 
an oblate intrinsic deformation with a triangular three-alpha-cluster configuration 
has been suggested. However, the HF+BCS calculation gives a spherical ground 
state. In the top-right portion of the figure we show the potential energy curve, 
which has only one minimum at the spherical shape. Other widely-used Skyrme 
forces of the SkM* and the SGlPJ also give the only minimum at sphericity. On 
the other hand, calculations using the Nilsson modeP^ and the Strutinski methocP^ 
give oblate ground states. An old Skyrme force SIl^'^' also gives an oblate minimum 
with 5 = —0.27 (the dash curve in the figure). 

Concerning ^'^C and ^'^Be, the experimental B(E2)| values indicate large defor- 
mations: |a2o|=0.82 for ^''C and 1.1 for ^''Be. However, the potential energy curves 
in Fig. 1101 have only the spherical minimum. The quantum fluctuation in shape may 
be able to account for the large B(E2)| since the curves are very soft toward prolate 
deformations. 

It is difficult to conclude definitely only from Fig. what deformation ^^C 
should bear in the HF approximation, since the optimal shapes of light nuclei are 
apt to be changed when effects beyond mean-field approximations are taken into 
account. The parity projection has been reported to be especially importanlP^ 
because the triangular three-alpha-cluster configuration violates the symmetry. 

Let us briefly report on the three-alpha-cluster linear-chain state of ^^C. By 
extending the potential energy curve in Fig. to larger 5^ we have found the first 
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excited minimum at 5 ~ 1.0, in good agreement with the result of the Nilsson model 
for the linear-chain state (5=1.1)!^ The excitation energy from our calculation is 
21 MeV. Though it is much larger than the experimental value of 7.654 MeV, the 
overestimation will be improved by the angular momentum projection. 




20 40 60 80 100 120 140 160 1 80 N 



Fig. 11. Same as in Fig.|Hl but for the hexadecapole deformation parameter 040. The grid indicates 
the locations of the magic numbers for spherical shapes. The results of FRDM is taken from 
RefP 

In Fig. 1111 we compare the axially symmetric hexadecapole deformation pa- 
rameter 040 between the present HF-I-BCS calculation and the FRDM result.^ A 
common feature between the results of the two models is that the sign of a^Q is 
positive in the bottom-left quarter of the major shells and negative in the top-right 
quarter. This behavior is perfectly understandable in terms of the density distribu- 
tion of pure-j single-particle wave functions. A systematic difference between two 
calculations is found in Z, < 50 where 1 040 1 is very small with exception of only 
several nuclei in the HF-I-BCS while it has certain sizes in many nuclei according to 
the FRDM. In particular, most of the light nuclei with N,Z < 20 have conspicuously 
large |a4o| in FRDM. For heavier nuclei, too, the magnitude of 040 is generally smaller 
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in HF+BCS than in FRDM. Further investigation is necessary to see whether these 
differences are originated in the different definitions of the shape parameters: In the 
HF+BCS they are calculated from the moments, while in the FRDM they are the 
input parameters to specify the shape of the single-particle potential. 

One of the advantages of mean- field methods over shell-correction schemes is that 
the protons and the neutrons do not have to possess the same radius and deformation. 
Indeed, it is one of the reasons to adopt the Cartesian-mesh representation rather 
than the expansion in harmonic-oscillator bases. In order to make the best use of 
this advantage, we have calculated the liquid-drop shape parameters separately for 
protons and neutrons for 1029 ground and 758 first-excited solutions. 




-0.3 0.3 0.6 -0.1 -0.05 0.05 0.1 

a,n OF NEUTRONS a.^ OF NFUTRONS 



Fig. 12. Difference of the shape parameters between protons and neutrons for the ground states 
of 1029 even-even nuclei (denoted by plus symbols) and the first-excited solutions (denoted by 
cross symbols) of 758 even-even nuclei for HF-(-BCS with Sill force. The left and the right 
portions show respectively the quadrupole deformation parameter 020 and the hexadecapole 
deformation parameter 040- 

The left-hand portion of Fig. El presents a comparison of the values of a2o 
between protons and neutrons. The r.m.s. value of the difference |fl2o"~^20°l 0.012, 
while the maximum difference occurs in the ground state of ^^Beio (^20*^ ~ 0.36, 
020° ~ 0.52), which is indicated by letter A in the figure. Symbols indicated by 
letters from B to F are the ground states of ^4Bei2 (020*^ = 0.44, 03™ = 0.59), 
fgAri2 (0.28, 0.18), f^Mgig (0.18, 0.27), i2Cri2 (0.36, 0.28), and fSAno (0.18, 0.11), 
respectively. Among these nuclei, ^'^Ar, ^^Cr, and ^^Ar are outside the proton drip 
line. The absolute value of the difference is smaller than 0.05 for 98.9 % of the 
solutions. 

The right-hand portion of Fig. is similar to the left-hand portion but for 
040. The r.m.s. value of the difference I040" — 040°] is 0.0063, while the maximum 
difference is found in ^4Bei2 (049" = —0.127, 04™ = —0.031), which is indicated with 
letter a. Symbols indicated by letters from b to g are the ground state solutions of 
^|Beio (a^o" = 0.065, 0^0° = 0.002), }|Mg6 (-0.004, -0.065), f^Mgio (0.063, 0.004), 
f^Nei2 (0.005, 0.063), f^Se (0.000, -0.045), and i|Ci2 (-0.048, -0.005), respectively 
Among these nuclei, '^^Mg and ^^S are outside the proton drip line. 
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To summarize, according to the results of our calculations with the Skyrme 
SIII force, the largest differences between the shapes of proton and neutron density 
distributions occur in the lightest nuclei. Except for these nuclei, the differences are 
not remarkably large. 




6 ti_ 



-0,3 0.3 0.6 -0.3 0.3 0.5 -0.1 0.1 

^20 (GROUND STATES) (GROUND STATES) 

Fig. 13. Correlation between quadrupole and hexadecapole deformation parameters. The left- 
hand portion shows the relation between a2o and 040 of each HF+BCS solution. Open circles 
are for ground state solutions while solid triangles are for the first excited solutions. The middle 
portion shows the relation between 020 of the ground state and 020 of the first excited state of 
a nucleus. An open circle corresponds to a nucleus. Only a kind of symbol is used. Clustering 
open circles may look like solid circles. The right portion is the same as the middle portion 
except it is for 040. 

In Fig. 1131 we investigate the correlations between deformation parameters. The 
left-hand portion of the figure shows that the ground states (open circles) are mostly 
with 0.2 < < 0.3 while the excited states (solid triangles) are with —0.3 < 020 < 
0. For prolate (a2o > 0) states 040 ranges from —0.1 to 0.1 while for oblate (020 < 0) 
states 040 covers a smaller interval from —0.05 to 0.07. For nuclei with 020 ~ 0, 040 
is also close to zero. 

The middle portion of Fig. [T^ shows that 020 of the ground state of a nucleus and 
020 of the first excited state of the same nucleus have a strong negative correlation. 
The most frequent case is that the ground state is prolate and the first excited state 
is oblate. There are opposite cases, too, but the number is much less. There are also 
many nuclei in which either ground or the first excited states is spherical (020 = 0). 

The right-hand portion of Fig. shows that 040 have more complicated corre- 
lation pattern between ground and the first excited states of a nucleus. 

A merit of extensive calculations is that it allows one to see global trends of 
nuclear properties over the entire area of the nuclear chart. As another example 
of the analysis of such a global trend, we have investigated the systematics of the 
energy difference between the oblate and the prolate solutions. In the top portion of 
Fig. 1141 the energy difference is plotted versus the neutron number. One can see an 
evident difference between below and above the = 50 shell magic: For nuclei with 
A^ < 50, the oblate solutions often have lower energies than prolate ones. For nuclei 
with A^ > 50, oblate ground states are very rare and found only in nuclei whose A^ 
is slightly smaller than a magic number. In the bottom portion, the same data is 
plotted versus the proton number Z. One can see that the situation is exactly the 
same for proton's shell effect. 

This abrupt change at A^ = 50 and Z = 50 seems to suggest that the dominance 
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Fig. 14. Energy difference between prolate and oblate solutions. Top and bottom portions show 
the same results for different abscissae. In the top portion, the abscissa is the neutron number 
and isotopes are connected with lines. It is possible to road off the proton number of each 
isotope chain from the direction of the "hand" put in each circle. See the legends in the top-left 
corner of each portion. In the bottom portion, the abscissa is the proton number and isotones 
are connected with lines. 



of prolate deformations in rare-earth and actinide nuclei may be related with the na- 
ture of the Mayer-Jensen major shell while the harmonic-oscillator major shell leads 
to an equal situation between oblate and prolate deformations. The Mayer-Jensen 
major shell is composed of the normal parity orbitals which have the same oscilla- 
tor quantum number and a unique parity orbital which is pushed down from the 
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next oscillator shell by the spin-orbit potential due to its large angular momentum. 
We have sugg 

estecPJ 

that one may ascribe the prolate dominance to the spin-orbit 

potential. 

Recently, we have performed Nilsson-Strutinsky calculations to answer this ques- 
tion.l^ The Nilsson potential contains a P potential, which simulates a square-well 
like radial profile of the mean potential, and a spin-orbit potential. The l'^ potential 
can reproduce the prolate dominance without the help of the spin-orbit potential, 
while the spin-orbit potential alone cannot reproduce the prolate dominance. How- 
ever, this does not mean the weakness of the effect of the spin-orbit potential. On the 
contrary, it interferes so strongly with the P potential that the ratio of the number of 
prolate nuclei to that of oblate ones oscillates with a large amplitude, making a pro- 
late dominance with the standard strength, reverting the situation to a oblate-favor 
one with half the standard strength, and again giving rise to a prolate dominance 
with vanishing strength. 

4.3. Radius 

The left-hand portion of Fig. ^1 presents a comparison of the liquid-drop radius 
of protons Rq^° with that of neutrons Rq'^^. 
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Fig. 15. The left-hand portion shows a comparison of the liquid-drop radius Ro of the proton 
distribution and that of the neutron distribution. The right-hand portion shows the definitions 
employed in this paper of the skin thickness and the halo radius. 



The r.m.s. value of the difference 



Tjneu 



Rq^°\ is 0.18 fm. Symbols indicated by 



letters from A to F are the solutions of f^Se (i^g^" = 3.3 fm, R^'" = 4.6 fm), fOSig 



(3.3, 4.2), i^Bei2 (3.9, 3.0), ^Cs (2.5, 3.4), 



22Ci6 (4.2, 3.3), 'and i|Beio (3.8, 2.9), 



respectively. Among the nuclei mentioned above, ^"Si is for the excited state while 
the others are for the ground states. ^^S, ^'^Si, and are outside the proton drip 
line. Hence the maximum difference for nuclei inside the drip lines occurs in the 
ground state of ^^Be. Incidentally, when we plotted the r.m.s. radii r^rns instead of 
Rq, we found that the resulting graph looked very similar to this figure except for a 
constant factor, Rq ~ (f)^''^ ?"rms- 

Concerning the difference of radial density distribution between protons and 
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neutrons, one should pay attention not only to volume properties represented by 
liquid-drop or r.m.s. radius but also to surface properties, which are sensitive to the 
difference at much lower densities and manifest themselves, e.g., in the thicknesses 
of nucleon skins and the radii of nucleon halos. The right-hand portion of Fig. 1151 is 
an illustration to explain our definitions of these two quantities. We regard that a 
point is in the proton skin if pp > 3/9n and Pp + Pn > 0.0016fm~^, where pp and pn 
are respectively the proton density and the neutron density at that point. A similar 
definition has been proposed in Ref.*^ The halo radius is defined as the radius at 
which angle-averaged mass density is 1.6 x 10~^ fm~^. 




-0.3-0.2-0.1 0.1 0.2 0.3 -0.3-0.2-0.1 0.1 0.2 0.3 



D/A D/A 

Fig. 16. Fermi levels as functions of the distance from the /3-stability hne D divided by the mass 
number A. The left and the middle portions are for protons and neutrons, respectively. The 
right portion explains the definition of D. 



As a preparation before discussing about skins and halos, let us first look at 
the relation between the Fermi levels and the distance in the (A^, Z) plane from the 
/3-stability line. As illustrated in the right-hand portion of Fig. 11(51 we choose to 
measure the distance D along a line of constant mass number. We assign a negative 
sign to D for nuclei with Z > N. As the definition of the /3-stability line, we employ 
a classical empirical formula, 

0.4^2 

N-Z = - , A = N + Z. (4-3) 

A + 200 ^ ^ 

The left- and right-hand portions of Fig. show respectively the Fermi level of the 
protons Ap and that of the neutrons An. Both ground-state and first excited solutions 
are plotted. The abscissa is the distance D divided by the mass number A. One can 
see a strong correlation in both portions. (Incidentally, the dots spread over larger 
area if the the abscissa is changed to DA~'^^^ or DA~^^^.) These figures suggest that 
one can measure the unstableness of a nucleus in two roughly equivalent ways. One 
may ask how high one of the Fermi levels is instead of asking how far the nucleus is 
located from the /3-stability line in the nuclear chart. 

Now, by using the next two figures, we will show that the skin grows monoto- 
nously and regularly as nucleons are added to the nucleus while the halo radius grows 
very slowly except near the drip lines, where it changes the behavior completely and 
expands very rapidly. 
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Fig. 17. Skin thickness. In the left portion the thickness of the proton (neutron) skin of the ground 
states is shown with the radius of open (solid) circles. Nuclei without either type of the skins 
are shown with the smallest dots. The staircase like lines represent two-proton and two-neutron 
drip lines. The right top portion shows the thickness of proton skin as a function of proton's 
Fermi level measured from the top of the Coulomb barrier as defined in the text. The right 
bottom portion shows the thickness of neutron skin versus neutron's Fermi level. 

The left-hand portion of Fig. El shows the thickness of the proton and neutron 
skins in the {N, Z) plane. Note that the absolute magnitudes of this quantity depend 
on our specific choice of the definition of the skins but the relative differences between 
adjacent nuclei are rather insensitive to the definition. The increase of skin thickness 
along isotope and isotone chains looks monotonous. If one would fold the graph in 
a line N = Z, one could easily find a general rule that the proton skin of a proton 
rich nucleus is larger than the neutron skin of its mirror nucleus. This fact can be 
ascribed to the Coulomb repulsion between protons. However, since the proton drip 
line is closer to the N = Z line than the neutron drip line is, proton skins are much 
rare than neutron skins if nuclei outside the drip lines are excluded. 

The top-right and the bottom-right portions of Fig.^lprssent the skin thickness 
of protons and neutrons, respectively. The Fermi levels are used as the abscissae. 
These figures include the data of not only those nuclei shown in the figure in the 
left-hand portion but all the calculated 1029 even-even nuclei. They also include not 
only the ground-state solutions but also first-excited solutions for 758 nuclei. Most 
of the dots are on the bottom line. From the figure in the right-bottom portion, 
one can see a general rule that the neutron's skin (by our definition) appears for 
An > — 8MeV and increases linearly with An at a rate of ^ 0.3fm/MeV. This rule is 
rather insensitive to the mass and applies both to ground and excited states. 

For protons, the Fermi level should be measured from the top of the Coulomb 
barrier. Otherwise the dots do not cluster in a narrow area. We estimate the height 
of the barrier top as 

1.44Z 

= 1.2^1/3 + Z^i? ^i? = 1.0fm, (44) 

where AR is a parameter to designate the shift of the location of the barrier top from 
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the liquid-drop surface. We assume a constant value of AR for all the nuclei. Small 
changes of AR only translate uniformly the plotted dots. We choose AR = 1.0 fm so 
that the horizontal location where the dots depart form the bottom line is roughly 
the same as in the graph for neutrons. The resulting behavior of the skin thickness 
of protons is very similar to that of the neutrons. 
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Fig. 18. The radius of the neutron (proton) halo plotted versus the thickness of the neutron 
(proton) skin (in the left (right) portion). Element names are printed for nuclei with large halo 
radius. An asterisk symbol put to the right shoulder of the element name indicates that the 
state is the first excited solution of the nucleus. 



Fig. 1181 shows "the thickness of the halo", i.e., the radius of the halo subtracted 
with the liquid-drop radius 1.2j4~^/^ fm. The abscissa is the thickness of neutron 
(proton) skin for the neutron (proton) halo. One can see that halos grow only slowly 
for skin thickness less than 1 fm but it expands rapidly after 1 fm. This sudden 
change is widely understood to be ascribed to the last occupied orbital's spatial 
extension due to its small binding energy. 
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Appendix A 

Cartesian Mesh and the Fourier Basis 

In this Appendix we explain the origin of the surprising high accuracy of nuclear 
mean-field calculations with apparently coarse Cartesian meshes. It is sufficient to 
consider a one-dimensional space because the essence is there. We put N mesh 
points at Xj = {i — (1 < i < N) in an interval < a; < Na as in the following 
illustration. 

a a a a ■ ■ ■ a 



Xi X2 X3 X4 



In the mesh representation, any function ^(x) defined in this interval is expressed in 
terms of a set of its values at the mesh points ■i/'i = "^(xi)- The equation to determine 
{ipi} is usually derived as a discrete approximation to the Schrodinger equation. An 
alternative point of view was presented by Baye and Heenen.EHJ They introduced a 
set of orthogonal basis functions fi{x), 



fi[x)=f{ , f{^) = — — , ko = —, l<i<N. (A-1) 

a J Asmafco^ A a 



The shape of the above basis function is shown in Fig. 1191 When ip{x) can be 




y=f«(x) 



mesh spacing: a 
box size ; 1 9a 
mesh points: ■> 



6 5a Wa T5a 

X 

Fig. 19. The graph of a Fourier basis function y — fsix) corresponding to the 8th mesh point for 
mesh spacing a and interval 19a. All the mesh points are indicated with open circles on the 
y = line. 

expressed as a linear combination of {/i(x)}, ipi is nothing but the coefficient of 

N 

i;ix)=^i;iMx), (A-2) 

i=l 

owing to a property of the basis, fi{xj) = 6ij. From this point of view, one can 
uniquely derive the equation to determine {ipi} from the variational principle in a 
space spanned by the trial wave functions of form HA-2() and can avoid the arbitrari- 
ness to choose a discrete approximation formula. 

The expression ()A-2|) is equivalent to a truncated Fourier expansion with a 
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boundary condition ip{x + Na) 



-1) 

N 



N+1 



'ip{x) because one can write 



(A-3) 



n=l 



This is the very reason why the Cartesian mesh representation provides high accuracy 
in nuclear calculations: The Cartesian mesh is associated with the Fourier or plane- 
wave basis, which is especially suitable to the atomic nucleus owing to its saturation 
property. 

Concerning integrals, the mid-point rule is exact for the product of two wave 
functions. 



Na 



ip*{x)(p{x)dx = ay~^ ip* 



(A-4) 



where ^p{x) and (j){x) are functions of form (|A-2|) and 0j = (j){xi). For other inte- 
grands like 'ip*{x)V{x)4){x), the mid-point rule similar to Eq. (|A-4p is not exact any 
more. Nevertheless, one should yet use the mid-point formula rather than high-order 
formulae like Eq. HC-2|) for the sake of high precision. 

On the other hand, the expressions for derivatives are not so simple. The ex- 
pression of the z^th derivative of the function defined by Eq. I)A-2|) at a mesh point 



Xi involves a full N x N matrix D, 



1 ^ 



(A-5) 



where the matrices for = 1 and 2 are expressed respectively as 





(1) 



-1) 



TT 



A^sin ^]Y^7r 



(i = j) 



(A-6) 



and 
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-IT 



-(-1) 



2n' 



N"^ tan ^]\rvr sin ^j^tt 



(i = j) 



(A-7) 



Since the evaluation of this expression seems to require a rather long computation 
time for the resulting precision, the code EV^ employs a finite-point approxima- 
tion formula like Eq. (|B-1|) to approximate Eq. ()A-5|) rather than sticking to the 
variational picture as completely as possible. 

The extension to the three-dimensional case is straightforward. A wave function 
is expressed by its values ipijk at mesh points {xi,yj, z^) = ((z — |)a, (j — i)a, (fc — i)a) 
as, 

'ip{x,y,z) = ^i^ijkfi{x)fj{y)fk{z). (A-8) 

ijk 
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It is worth noting that, for polar and cylindrical coordinates, the mesh points 
become non-uniform in order to make the associated bases the eigenfunctions of the 
kinetic energy. If one takes the radial mesh points equidistantly, the spacing must be 
finer than that of the three-dimensional Cartesian mesh. Meshes with 1 fm spacing 
are too coarse for the radius grid but fine enough for three-dimensional Cartesian 
mesh. 



Appendix B 

Finite-Point Formulae for Numerical Differentiation 



The vth derivative of a function f{x) can be evaluated approximately by the 
values of the function at n mesh points as, 



1 



{n-l)/2 

E 

:={n-l)/2 



» 



c„//(x + ai) + 0(a- 



(B-1) 



where n > 3 is an odd number while [j]e means the maximum even number not 
exceeding j. The coefficients c^^ and bn^ can be determined by approximating f{x) 
with the polynomial in x of degree n — 1 which coincides with f{x) at the n points. 

Coefficients up to the six-point formula are given in table 25.2 of Ref.^l We 
have utilized a symbolic computation software mathematic^^ to obtain the values 
of the coefficients of higher-order formulae, which are given in Table 

We show in Fig. 20] the results of our calculations with which we examine the 
precision of the kinetic energy calculated using Eq. (in^l) . The results are for one- 
dimensional case but can be easily applied to three-dimensional cases because the 
kinetic energy is divided into x, y and z components. 
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Fig. 20. Precision of finite-point numerical differentiation formulae. 



In the left-hand portion of Fig. I2U1 the ratio of the numerically calculated value 
of Jq^^'' sin kx-^ sin kxdx to its exact value —k-K is plotted as a function of the 
wave number k. The mesh spacing a (expressed as Ax in the figure) is taken as 1 
fm. We change k from to 2tt for an aesthetic reason, though k > vr/a = vr fm~^ 
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Table I. Coefficients of the n-point approximation formulae for the i/th derivative given in Eq. l|B-lfl . 
^^ni foi' negative i can be obtained using a relation cJ^'^L^ = {~'^Y'^ni ■ '^no ~ ^ ^^o'' ^'^'^ ^■ 
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19 


12252240 


11027016, 4009824, 1559376, 539784, 154224, 34272, 5508, 567, 28 


21 


232792560 


211629600, 79361100, 32558400, 12209400, 3907008, 1017450, 205200, 
29925, 2800, 126 



Second derivative {u = 2) 



n 




lcjf/|;i = 0,l,2,... ,i(n-l) 


3 


1 


2, 1 


5 


12 


30, 16, 1 


7 


180 


490, 270, 27, 2 


9 


5040 


14350, 8064, 1008, 128, 9 


11 


25200 


73766, 42000, 6000, 1000, 125, 8 


13 


831600 


2480478, 1425600, 222750, 44000, 7425, 864, 50 


15 


75675600 


228812298, 132432300, 22072050, 4904900, 1003275, 160524, 17150, 900 


17 


302702400 


924708642, 538137600, 94174080, 22830080, 5350800, 1053696, 156800, 
15360, 735 


19 


15437822400 


47541321542, 27788080320, 5052378240, 1309875840, 340063920, 
77728896, 14394240, 1982880, 178605, 7840 


21 


293318625600 


909151481810, 533306592000, 99994986000, 27349056000, 7691922000, 
1969132032, 427329000, 73872000, 9426375, 784000, 31752 



is meaningless from the point of view of the Fourier expansion. In the numerical 
calculation, the second derivative of sinkx is approximated by that of the (2n + 1)- 
point Lagrange polynomial, which has degree n and coincides with sm kx at x, x i 
o, • • • , X lb na. The number of points is chosen as n = 3, 5, • • • , 25. The integral, too, 
is evaluated numerically but with much higher precision than that of the numerical 
differentiation. Hence the plotted quantity is the error of numerical differentiation, 
not that of integration. From this result one can see that the value of kinetic energy 
is always underestimated and the error becomes smaller by employing higher-order 
formulae. 

The middle portion of Fig.IJOJis the same as the left-hand portion except that the 
ordinate is the relative error | (ratio) — 1| in logarithmic scale. One can read off from 
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the figure, for example, tliat the relative error of the kinetic energy at k = 1.4 fm~ 
is 3.0 X 10"'^ for the 9-point formula of the second derivative. The 9-point formula is 
used in the calculations presented in this paper. This error roughly accounts for the 
error of the total energy of ^^'^Pu, which is 0.5%. One can conclude that the main 
source of the error of the total energy is the kinetic energy. 

In the right-hand portion of Fig. I20| the error of numerically evaluated values 
of the second derivative of a function at x = ^ is shown as a function of 

the mesh spacing a (called Ax in the figure). The numerical evaluation means the 
usage of the (2n + l)-point Lagrange polynomial which coincides with the function 
at x = i, i lb a, • • • , | ±na. The solid curve labeled as "Fourier" indicates the 
error of the Fourier interpolation, which coincides with the original function e'^^^"^ 
at infinite number of discrete points of x = ^ + na (n = 0, ±1, ±2, • • • ). For plane 
waves with k < tt/o, the Fourier interpolation gives the exact value. For functions 
involving higher wave-number components, as the present Gaussian function, the 
Fourier method does not result in the exact value. Nevertheless, as one can see from 
the figure, it is still much more precise than polynomial-based formulae. 

In order to evaluate the derivative values at points within ^{n — l)a from the 
boundaries, one needs the values of the function beyond the boundaries. In order to 
obtain them one has to assume either vanishing (reflection antisymmetric) or periodic 
boundary conditions. A much easier method is to assume the function to be zero 
beyond the boundaries. However it is justified only when the function is well-damped 
before reaching the boundaries. Otherwise it produces a discontinuity of derivative 
values at the boundaries and invalidate the n-point formulae for n > 3. Usage of a 
cavity whose shape is not a box (rectangular parallelepiped) can be justified only for 
functions well-localized around the center of the cavity. 

Appendix C 

Finite-Point Formulae for Numerical Integration 

The integral of a function f{x) over an interval < a; < L can be evaluated by 
the values of the function at n -|- 1 mesh points, Xi = ai {0 < i < n, a = L/n), hy 
approximating the function with a polynomial in x of degree n coinciding with f{x) 
at the mesh points. Such formulae are called (n-|- l)-point Newton-Cotes closed type 
integration formula and are expressed as 

fn.)4.^a±^nai) + Oia^*>). *={^ (--j;) . (CI) 

We have calculated the values of the coefficients a^i and /?„ in a similar method as 
we calculated those of the differentiation formulae. The results are given in Table ITII 
for some odd values of n + 1. 

Two-, three-, and four-point formulae are called respectively the trapezoidal 
rule, Simpson's ^ rule, and Simpson's | rule. The errors of 2/c-point and {2k — 1)- 
point formulae are of the same order in a, ©(a^^^+i), for 2k > 4. For this reason we 
listed in Table ITll onlv 2- and odd-point formulae. The coefficients of 4- to 11-point 
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Table II. Coefficients of Newton-Cotes closed type (n + l)-point integration formulae. Values for 
i > |(n + 1) can be obtained from a relation a^i ~ an.n-i- 



n+1 


/3n 


a„i; i = 0, 1, • • • ,^n+l 


3 


3 


1, 4 


5 


45 


14, 64, 24 


7 


140 


41, 216, 27, 272 


9 


14175 


3956, 23552, -3712, 41984, -18160 


11 


299376 


80335, 531500, -242625, 1362000, -1302750, 2136840 


13 


5255250 


1364651, 9903168, -7587864, 35725120, -51491295, 87516288, -87797136 


15 


2501928000 


631693279, 4976908048, -5395044599, 24510099488, 
-46375653541, 88410851312, -117615892611, 136741069248 



formulae agree with Eqs. 25.4.13 - 20 of Ref.l^ There is no difficulty to obtain the 
coefficients of formulae which involve more number of points. It should be kept in 
mind, however, that formulae for n -|- 1 >= 9 (n + 1 = 10 is an exception) include 
negative coefficients, which cause a loss of numerical precision due to cancellations 
between the contributions of negative and positive coefficient points. This problem 
becomes more serious for larger n. 

In order to apply the (n -|- l)-point formula to a long interval, one divides the 
interval into m segments and applies the (n-l-l)-point formula to each of the segments, 
i.e., 

fmna m—1 C n f \ 

- E E f / (c^ + + - ? :r„? ■ 

(C-2) 

Note that points x = jna {1 < j < m — \) appear twice in the double summation. 
The reason for 5' = 6 — 1 \s that the error of the integral over the interval length 
L = mna is m = L/na oc times as large as the error on the right-hand side on 

Eq. jnu). 

We have calculated the precision of these formulae as a function of the mesh 
spacing a for three types of integrals, 

/ -T, — dx, / dx, / e ^ dx, C-3 

Jo 2x Jq 1 + x^ Jq 

as examples of three different boundary situations at x = 27r. The obtained (absolute 
value of the) relative errors of the integrals are shown in logarithmic scale versus the 
mesh spacing a (= Ax) in Fig. 1211 The calculations are defined only for discrete 
points at a = 27r/mn. The lines drawn in the figure connect such discrete points. 

In the left-hand portion of Fig. 1211 the error of the n-point formula behaves 
as expected from Eq. HC-2|) for every value of n. More-point formulae give quicker 
decreases of the error as the mesh spacing decreases. The noise-like error of magni- 
tude 10-^3 to 10"^^ originates in the round off of real numbers into double precision 
(16-digit) floating numbers. 

However, this is not the only situation. We now show that the convergence to the 
exact value is much faster for some kind of integrands and, what is more interesting is. 
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that such fast convergences become more outstanding by using low-order few-point 
formulae than high-order many-point ones. 

The right-hand portion of Fig. shows the results for a different integrand 

2 

. In this case, one can see that fewer-point formulae give more accurate results 
than more-point formulae. The convergences are generally much faster than those 
in the left-hand portion of the figure. 

The origin of this unexpected situation apparently in contradiction with Eq. HC-2() 

2 

is the approximate periodicity of the integrand. The function in an inter- 

val [0, 27r] can be regarded as a half wave length portion of a periodic function, 
which is constructed using this portion of the function by repeating reflections at 
the boundaries x = and 2tt. At x = the function is reflection symmetric orig- 
inally. For X > 2tt the function is redefined by a reflection symmetry requirement 
f{2TT + x) = f{2TT — x). This reflection makes practically no discontinuities at x = 2tt 

2 

(or between x = ±27r because the period is 47r) because and its derivatives (of 
not extremely high order; the graph includes up to 11-point formula and thus only 
up to 11th derivative concerns) are close to zero at x = 2Tr. 

The integral over [0, 27r] is a half of the integral of thus constructed smooth pe- 
riodic function over the wave length An. It is known that the best numerical formula 
for integrals of periodic functions over a wave length is the mid-point rule. The rea- 
son is often explained in terms of an asymptotic expansion. A more quantitatively 
accurate explanation can be given in terms of Fourier expansions, as explained in 
Appendix El 

Periodic functions composed of wave numbers less than vr/a can be exactly in- 

2 

tegrated with the mid-point rule. The function is an approximately periodic 
function and its large-wave-number components are non-zero but very small. There- 
fore the mid-point rule is the most precise formula for this integrand. The 2-point 
formula is the same as the mid-point rule in the present situation 

The (n -|- l)-point formulae can be regarded as a summation of n similar results 
of 2-point formula with mesh spacing of na. The curve of the error of the 3-point 
formula resembles to that of the 2-point formula except it is shifted leftward by a 
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factor (2 — l)/(3 — 1) = 2- Similarly, the horizontal distance between the curves for 
3- and 5-point formulae is also calculated as (3 — l)/(5 — 1) = |. 

We have searched for an integrand which gives rise to an intermediate situation 
between the left-hand and the right-hand portions of Fig. 1211 A function which 
serves for this aim is {l - (x/27r)2} / (1 + The middle portion of Fig. HU shows 
the error of the integral for this integrand. One can see that the situation is like that 
of the left-hand (right-hand) portion of the figure for small (large) a (= Ax). The 
most precise formula changes with the mesh spacing: It is the 2-point formula for 
a > 0.3, the 3-point formula for 0.08 < a < 0.3, and the 5-point formula for a < 0.08. 

This change can be understood as follows. With large a, the precision is so low 
that the discontinuity at x = 27r looks too small to invalidate the approximation of 
the periodicity. Consequently, the error behaves like a parabola as in the right-hand 
portion of the figure. With small a, the precision is so high that the discontinuity 
looks large enough to violate the periodicity. Thus the error behaves linearly as in 
the left-hand portion of the figure. 

In the nuclear mean-field calculations, the integrands are not exactly periodic. 
However, the situation resembles that of the right-hand portion of Fig. |^ for nuclei 
whose Fermi levels are not very shallow, because the wave functions are well damped 
before reaching to the boundaries. This is why we employ the mid-point rule in our 
calculations. 

There are also modified types of Newton-Cotes formulae, in which the coefficients 
are different from 1 only near the boundaries: 

I mdx = aY^^'^fia^) + Oia^-^'), 5' = [l ^^^^^^ , (04) 

with 

Oi'N,n,i =Ol'N,n,N-i (0<i<n), (C-5) 

a'Nni = P'n {n + l<i<N-n-l). (C-6) 

We temporarily call the above equations the boundary-correction type formulae in 
this paper. We have calculated the values of the coefficients a^„j and P'^ in a similar 
manner as we treated Eq. ()CT|) . Table HTll gives the values of the coefficients. For 
(ra + l)-point formula, the coefficients are different from one only for n + 1 points near 
each ends of the integration interval. Only in the 2-point formula, which is identical 
to the extended trapezoidal rule, the end points alone have a non-unit coefficient. 
We give both 2A;-point and (2A; — l)-point formulae in the table. Although their errors 
are of the same order in a, the former seems slightly more precise than the latter in 
some test calculations. We have found only the 3-point formula in literature. It is 
Eq. (4.1.14) of Ref.ESJ Although the formulae are identical, the derivation is quite 
different. 

Boundary-correction type formulae will be useful to treat deeply bound, weakly 
bound, and positive-energy orbitals in a single framework. With the formulae, the 
deeply bound orbitals are treated with the high precision of the mid-point rule for 
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Table III. Coefficients of Boundary-correction type (n + l)-point integration formulae. 



7Z -r 1 




CtNni', * — 0, 1, • ■ • , n 


2 


2 


1, 2 


3 


24 


9, 28, 23 


A 
4 




o, ol, zU, zo 


5 


1440 


475, 1902, 1104, 1586, 1413 


6 


1440 


459, 1982, 944, 1746, 1333, 1456 


/ 




oD/yy, i(DD4o, o4ooi, iMy84, 5y4o(, iouyoo, iiyooo 


8 


120960 


35584, 185153, 29336, 220509, 46912, 156451, 111080, 122175 


9 


7257600 


2082753, 11532470, 261166, 16263486, -1020160, 12489922, 5095890, 
7783754, 7200319 


iU 


/ zo / duu 


on'iAfiOK iio^ip;^ioo I/171/I/10 onQnt^oQQ vaq/ioqq ^Q^^A^\^f\ iriPiQiQQ 
zUo4dzO, llyDODzz, -14/144Z, zUoUDzoo, -^Uo4zoo, io004UoU, lUOoioo, 

9516362, 6767167, 7305728 


11 


958003200 


262747265, 1637546484, -454944189, 3373884696, -2145575886, 
3897945600, -1065220914, 1942518504, 636547389, 1021256716, 


12 


958003200 


257696640, 1693103359, -732728564, 4207237821, -3812282136, 
6231334350, -3398609664, 3609224754, -196805736, 1299041091, 
896771060, 963053825 


13 


5230697472000 


1382741929621, 9535909891802, -5605325192308, 28323664941310, 
-39865015189975 53315913499588 -41078195154304 39099895874876 
-13155015007785, 12465244770050, 3283609164916, 5551687979302, 
5206230892907 


14 


5230697472000 


1360737653653, 9821965479386, -7321658717812, 34616887868158, 
-48598072507095, 81634716670404, -78837462715392, 76782233435964, 
-41474518178601, 28198302087170, -3009613761932, 7268021504806, 
4920175305323, 5252701747968 


15 


62768369664000 


16088129229375, 121233187986448, -109758975737401, 
502985565327936, -823993097730133, 1461175500619600, 
-1668277571373981, 1746664157478912, -1219701572786787, 
819644306759856, -276710928642475, 174692180291008, 
37176466501689, 66395850785776, 62528161418177 



periodic integrands while the weakly bound and positive-energy orbitals are treated 
with less precise but still high-precision (n -|- l)-point formula. 
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